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Stable embedded solitons are discovered in the generalized third-order nonlinear Schrodinger equa- 
tion. When this equation can be reduced to a perturbed complex modified KdV equation, we devel- 
oped a soliton perturbation theory which shows that a continuous family of sech-shaped embedded 
solitons exist and are nonlinearly stable. These analytical results are confirmed by our numerical 
simulations. These results establish that, contrary to previous beliefs, embedded solitons can be 
robust despite being in resonance with the linear spectrum. [PACS: 42.65.Tg, 05.45.Yv.] 



Pulse propagation in optical fibers attracts a lot of at- 
tention these days. Pico-second pulses are well described 
by the nonlinear Schrodinger (NLS) equation which ac- 
counts for the second-order dispersion and self-phase 
modulation. But for femtosecond pulses, other physi- 
cal effects such as the third-order dispersion and self- 
steepening become non- negligible. The physical model 
which incorporates these additional effects is [1,2] 

i [V'z + f^ltpl'tpT + fi{\tljf)ri^] + i^rr + \i^fi> + itpTTT = 0. 

(1) 

Here ip is the envelope of the electric field, z is the 
distance, r is the retarded time, P and are the self- 
steepening coefficients [1]. All quantities have been nor- 
malized. For optical pulses, (3 = jj, [1,2]. In this paper, 
we allow (3 and /i to be different for the sake of mathe- 
matical analysis. The Raman effect, which is dissipative 
in nature, is also non- negligible [1,2]. It is not included 
in the model (1) because we want to focus our attention 
on the other physical effects in this paper. 

The third-order dispersion term in Eq. (1) is signifi- 
cant because it qualitatively changes the linear dispersion 
relation of Eq. (1). Its effect on the NLS soliton is to 
generate continuous-wave radiation and causes the soli- 
ton to decay [1-5]. However, solitary waves of Eq. (1) 
which are embedded inside the linear spectrum do ex- 
ist in certain parameter regimes [6-9], and such waves 
are called embedded solitons [10,11]. To see why a soli- 
tary wave of Eq. (1) has to be an embedded soliton, we 
substitute solitary waves of the form ^'(r — Vz)e^'^^ into 
(1), where velocity V and frequency A are constants. We 
readily find that for any frequency A, the linear equa- 
tion for ^E* allows oscillatory solutions. Thus, all A lies in 
the continuous spectrum of Eq. (1), hence the solitary 
wave must be an embedded soliton. Stability of these 
embedded solitons is clearly an important issue. Previ- 
ous analytical studies have shown that if embedded soli- 
tons are isolated in a conservative system, they are at 
most semi-stable, i.e., the perturbed soliton persists or 
decays depending on whether the initial energy is higher 
or lower than that of the embedded soliton [10,12,13]. A 
physical explanation for it is as follows [10]. If the initial 
state has energy higher than the embedded soliton, it just 



sheds extra energy through tail radiation and asymptoti- 
cally approaches this embedded soliton; but if the initial 
state has lower energy than the embedded soliton, the 
energy loss (through radiation) drives the solution away 
from the embedded soliton, resulting in instability. How- 
ever, embedded solitons in Eq. (1) are not isolated: they 
exist as a continuous family, parameterized by their en- 
ergy [8,9]. Thus, the above physical argument for semi- 
stability does not apply. Numerical results suggest that 
the families of double-hump embedded solitons in Eq. 
(1) are still semi-stable [9], but the family of single- hump 
embedded solitons are fully stable [14]. 

In this paper, we show that single-hump embedded 
solitons in Eq. (1) are fully stable when [/? — 6] ^ 1 
and ]/i] <C 1. To our knowledge, this is the first time sta- 
ble embedded solitons are rigorously established in the 
literature. This result indicates that embedded solitons 
can be as robust as conventional (non-embedded) soli- 
tons, contrary to the previous belief. The method we 
will use is to develop a soliton perturbation theory for 
the (integrable) complex modified KdV (CMKdV) equa- 
tion, supplemented by numerical simulations. 

We first employ the variable transformation 
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(2) 



upon which Eq. (1) becomes 

Uz + uj^j + 6]u]^uj — ia\u\'^u — ^{\u\'^)^u, (3) 

where a = (6 — /3)//3, and 7 = 6/x//3. If a = 7 = 0, the 
above equation is the CMKdV equation which is com- 
pletely integrable by the inverse scattering method [15]. 
It admits sech-shaped soliton solutions whose amplitudes 
and velocities are free parameters (see below). In this pa- 
per, we consider the case \a\, \^\ <^ 1, i.e., \(3 — 6\ ^ 1 
and \fi\ <ti 1. In this limit, soliton evolution in the per- 
turbed CMKdV equation (3) can be studied by a soli- 
ton perturbation theory. For this purpose, we denote 
a — eao,j — ejo, where \e\ ^ 1. When e = 0, Eq. (3) 
has the following soliton solutions 



uoiz,0 = Uo{e)e 



iXz 



^-Vz, 



(4) 



where 



1 



Uoie) = r8ech.rde'^\ V = - 3A;^ A : 



-2fc(r^+A;^). (5) 



Here amplitude r and frequency k arc free parameters. 
When ^ |e| <C 1, these sohtons will deform due to per- 
turbations on the right hand side of Eq. (3), and their 
amplitudes and frequencies will undergo slow evolution 
with respect to z. Below, we use the soliton perturbation 
theory to derive this evolution. 

First, we write the solution in the form 



'Io^'^''U{z,e), = ^- f Vdz. 

Jo 



(6) 



Substituting this form into (3), we find the equation for 

Uiz,e) as 

U, + i\U - VUe + Ueee + Q\U\^U0 = eF{U), (7) 
where 

F{U)=[iao\U\^-M\U\MU- (8) 
Next, we expand the solution U into a perturbation series 

[/(z, 9) = Uo{e) + eUi{z, e) + €^U2{z, 6) + ..., (9) 

where Uq is given in (5). When this series is substituted 
into Eq. (7), at order 1, the equation is satisfied automat- 
ically. At order e, we obtain the linear inhomogeneous 
equation for Ui as 



F{Uo) 



(10) 



where C is the linearization operator of Eq. (7), and su- 
perscript "*" denotes the complex conjugation. At initial 
distance z = 0, Ui = 0; when z ^ 1, Ui approaches a 
steady state with continuous-wave tails at infinity. Gen- 
eration of continuous- wave tails is a distinctive feature of 
embedded solitons under perturbations [10]. These tails 
have the form /le"^**^^, where h is the tail amplitude. 
Due to the Sommerfeld radiation condition, these tails 
can only appear at 9 -C —1, not at ^ » 1. In other 
words, 



lim Ui 
+00, e<-i 



he' 



-2ik0 



lim Ui 
+00, e»i 



0. (11) 



One of the key steps in the soliton perturbation theory 
is to determine the amplitude h of continuous-wave tails 
in Ui. This can be done by solving Eq. (10) directly, 
starting with the zero initial condition [13]. The way to 
do it is to expand the solution into the complete set of 
eigenfunctions for the linearization operator £, which are 
just the squared eigenstates of the Zakharov-Shabat sys- 
tem with the soliton potential [16]. This we have done. 
But a much simpler way to derive h is to just consider 
the steady-state solution Ui . This allows us to drop the z 
derivative in Eq. (10), and use only solvability conditions 
to obtain h. This idea has been used before [12]. 



To pursue this latter approach, wc need the bounded 
eigenfunctions of the adjoint linearization operator 
with zero eigenvalue. Here the inner product used to 
define an adjoint operator is 



f{0fg{e)de, 



(12) 



where the superscript "T" represents the transpose of a 
vector. Under this inner product, the adjoint operator 
£^ can be readily obtained. Eigenfunctions of jC"^ are 
simply a different set of squared Zakharov-Shabat eigen- 
states [16]. With this in mind, we can easily show that 
has four bounded eigenfunctions for zero eigenvalue 
— two localized and two non-localized. The localized 
(discrete) eigenfunctions are 



$1 = 



TT* 
-Uo0 



$2 



US 



(13) 



They are associated with the translational and phase in- 
variances of solitons in the CMKdV equation. The non- 
localized (continuous) eigenfunctions are 
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r2sechV6'e-4*'^^ 
{rtanhrO + 2,ikf'e 



-2ike 



(14) 



(15) 



Now we take the inner product between Eq. (10) (with 
no z derivative) and each of the above four eigenfunc- 
tions. In view of the asymptotics (11), we see that the 
first two inner-product equations are simply 



F{Uo) 
F*{Uo) 



$2 



FiUo) 
F*{Uo) 



= 0, 



(16) 



which are satisfied automatically when the form (8) for 
F{U) is utilized. The last two inner-product equations 

arc equivalent. After simple calculations, these equations 
give a formula for the amplitude h as 



7 1 • / ^ ^ , Snk 

^ = -4^^("° + 27ofc)^:^sech— 



-TT-. (17) 



When the tail amplitude h has been derived, we can 
proceed to order e^. Again, by applying the solvability 
conditions for solution U2, the dynamical equations for 
soliton parameters r and k will be obtained (see [12] for 
an example). An alternative way to derive these dynam- 
ical equations, which is simpler and physically more in- 
sightful, is by using conservation laws [12]. In this latter 
way we will proceed. 

Eq. (3) admits the following two conservation laws: 



h 



f 

o — c 



(18) 
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J2 = y I iia(«{U* - uuj) - 7|u{ 1^ + 7(1 + i7)|u|''| d^. (19) 

When Eq. (6) and the perturbation series (9) are substi- 
tuted into these conservation laws and terms up to order 
retained, the first conservation law becomes 




{\Uof + €{UoU; + U;Ui) + e'^iUoU^ + !7*!72) + e^|(7i|^} = 0. 



(20) 

To proceed further, the following facts are noted: (i) |C/oP 
changes slowly (on the distance scale e^z); (ii) the Ui 
field is driven by the inhomogeneous term in Eq. (10) 
and quickly becomes stationary in the soliton {Uq) re- 
gion; the U2 field is similar; (iii) the Ui field develops 
a continuous-wave tail of amplitude h on the left-hand- 
side of the soliton [see Eq. (11)]; this tail travels at its 
group velocity Vg = — 12/c^; in the moving coordinate ^, 
its relative velocity is Vg — V ~ — (r^ + 9fc^) [see Eq. 
(5)]. The fact (ii) means that the z derivatives of the 
integrals of the second and third terms in Eq. (20) are 
zero and can be dropped. The fact (iii) means that the 
z derivative of the integral of the last term in (20) is 
e^{r^ + 9k^)\h\^. When these results and the Uq formula 
(5) are utilized, the first conservation law (20) gives the 
dynamical equation for the soliton amplitude r. Similar 
calculations for the second conservation law produces the 
dynamical equation for the soliton frequency k. When 
the h formula (17) is substituted, and (eao, 670) replaced 
by the physical parameters (a, 7), the dynamical equa- 
tions for r and k finally take the form 

fz = + 2^^)'^"' + 9A:^)sech^^, (21) 

^ = -^7r^(a + 2'yk)[3ak — 7(r^ -|- 3k^)]{r^ + 9fc^)sech^^^, 

(XZ 'jZii ^1 

(22) 

and they govern the soliton evolution under small per- 
turbations \a\, I7I ^ 1 in Eq. (3). These two equations 
are the main analytical results of this paper. 

Dynamical equations (21) and (22) admit a continuous 
family of fixed points: 

a 

k = — — , r is arbitrary. (23) 
27 

These fixed points correspond to a continuous family of 
embedded solitons (4) in the original wave equation (3). 
Their amplitudes r are arbitrary, while their speeds V de- 
pend on r according to Eq. (5). It is easy to check that 
these embedded solitons are precisely the ones reported 
in [8] . We have further determined the linear and nonlin- 
ear stabilities of these fixed points in equations (21) and 
(22), and found that they are all stable, both linearly 
and nonlinearly. This means that these sech-shaped em- 
bedded solitons in the perturbed CMKdV equation (3) 
are both linearly and nonlinearly stable, despite the fact 



that such solitons reside inside the continuous spectrum 
of the wave system and are in resonance with the contin- 
uous waves. They mark a distinct contrast with isolated 
embedded solitons and the continuous family of double- 
hump embedded solitons in other physical systems, which 
were shown to be semi-stable [9,10,12,13]. 

The simplest way to demonstrate the linear and non- 
linear stabilities of fixed points (23) is by drawing a phase 
portrait in the (r, k) plane. Such a phase portrait with 
a = —0.1 and 7 = 0.1 is shown in Fig. 1. At these param- 
eter values, the fixed points are the k = ^ line (dashed). 
Clearly, we see that these points are nonlinearly stable, 
and they attract a large region of initial parameters. At 
other a and 7 values where < 0, the phase portraits 
are qualitatively similar. If 0:7 > 0, the phase portrait is 
just the one of 0:7 < with k replaced by —k. 



k 




r 

FIG. 1. Phase portrait of the dynamical system (21) and 
(22) for physical parameters a = —0.1 and 7 = 0.1. 



Lastly, we confirm the above analytical results by di- 
rect numerical simulations of Eq. (3). Our numerical 
scheme is the pseudo-spectral method along the ^ direc- 
tion, and the fourth-order Runge-Kutta method along 
the z direction. Damping boundary conditions have also 
been used to filter out the energy radiation emitted into 
the far field. A number of numerical simulations of Eq. 
(3); starting with the soliton solution (4) at z = 0, 
have been performed for various small a and 7 values. 
Good agreement between the numerics and the analysis 
has been observed. In addition, as a and 7 approach 
zero, the difference between the leading-order perturba- 
tion results [Eqs. (21) and (22)] and numerical values 
approaches zero, meaning that our perturbation theory 
is asymptotically accurate. For illustration purpose, we 
fix the physical parameters a = —0.1 and 7 = 0.1, and 
choose two sets of initial amplitude and frequency values, 
(ro, A:o)=(1.5, 1) and (1.5, 0), which are on the opposite 
sides of the line of fixed points in Fig. 1 . The simulation 
results are displayed in Fig. 2. For (ro, /co)=(1.5, 1), the 
pulse evolution is plotted in Fig. 2(a). Due to perturba- 
tions in Eq. (3), the speed of the pulse slowly increases. 
Thus the pulse turns around and eventually approaches 
a steady positive speed. In Fig. 2(b), the field profiles 
at two distances are plotted. As predicted, a continuous- 
wave tail develops on the left of the embedded soliton. 
Note that the tail amplitude decreases with distance z, 
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meaning that the embedded soHton is stabiHzing [10]. A 
direct comparison between the numerics and the leading- 
order perturbation theory in the phase plane is shown in 
Fig. 2(d). We see that the numerical pulse frequency k 
indeed approaches the theoretical value ^. 

When {ro,ko)—{l-5, 0), the numerical pulse evolution 
is displayed in Fig. 2(c). In this case, the soliton slows 
down under perturbations, contrary to the previous case. 
Examination of the radiation field indicates that the 
continuous-wave tails emitted to the left of the soliton 
are also lessening, indicating that the embedded soliton 
is stabilizing again as the theory predicts. Comparison 
in Fig. 2(d) between the numerics and the theory for this 
case also shows good agreement. 

Physically, why are embedded solitons in Eq. (3) fully 
stable? There are two major reasons. First, these soli- 
tons exist as a continuous family. In other words, in the 
neighborhood of each embedded soliton, there are other 
embedded solitons nearby which have higher or lower en- 
ergy. Thus when perturbed, an embedded soliton may 
always relax into an adjacent one, similar to the NLS- 
solitons. The second reason is that these embedded soli- 
tons are single-humped. Recall that double-humped em- 
bedded solitons in Eq. (1) also exist as a continuous fam- 
ily, but they are not stable [9]. This is not surprising, as 
the instability of multi-hump solitons is well documented 
in the literature. Embedded solitons in Eq. (3), however, 
are single-humped. Thus they could be fully stable as we 
have shown in this paper. 

We note that for ultra-short pulses, (3 — fi; while in this 
paper, \a\, \^\ <^ 1, i.e., |/3 — 6| ^ 1 and ^ 1. In the 
physical case /3 = /i, a continuous family of single-hump 
embedded solitons also exists [8]. In addition, Gromov, et 
al's numerical computations show that a sech pulse tends 
to one or a few single-hump embedded solitons [14]. This 
numerical evidence, together with the above analytical 
results and physical arguments, strongly suggests that 
embedded solitons in this physical case are also stable. 
A rigorous proof of this full stability in this physical case 
can be provided by an elaborate internal-perturbation 
method [13], which we will pursue in the near future. It 
is also noted that for ultra-short pulses, the Raman effect 
is non-negligible. But since embedded solitons without 
the Raman effect are robust and stable, we can expect 
that the Raman effect on embedded solitons would be 
also a frequency downshift, similar to that on NLS soli- 
tons [1,2]. 

In summary, we have discovered fully-stable embed- 
ded solitons in a physical model relevant for ultra-short 
pulse propagation in optical fibers. This finding dispels 
previous skepticism about the observability of embedded 
solitons. We expect this work to have important impli- 



cations to ultra-short pulse propagation. It would also 
stimulate further search for stable embedded solitons in 
other physical systems. 

This work was supported in part by a NASA grant. 
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FIG. 2. Numerical evolution of initial solitons (4) with 
To — 1.5, fco = 1 [in (a, b)] and [in (c)] under perturbations 
of Q = —0.1 and 7 — 0.1 in Eq. (3). In (d), numerical and 
analytical results are compared in the phase plane for both 
cases of fco = 1 and 0; solid lines: analytical; dots: numerical. 
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